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Abstract-V^ describe  a  new  VNIR-SWIR  atmospheric 
correction  method  for  multi-  and  hyperspectral  imagery,  dubbed 
QUAC  (Quick  Atmospheric  Correction)  that  also  enables 
retrieval  of  the  wavelength-dependent  optical  depth  of  the  aerosol 
or  haze  and  molecular  absorbers.  It  determines  the  atmospheric 
compensation  parameters  directly  from  the  information 
contained  within  the  scene  using  the  observed  pixel  spectra.  The 
approach  is  based  on  the  empirical  finding  that  the  spectral 
standard  deviation  of  a  collection  of  diverse  material  spectra, 
such  as  the  endmember  spectra  in  a  scene,  is  essentially  spectrally 
fiat.  It  allows  the  retrieval  of  reasonably  accurate  reflectance 
spectra  even  when  the  sensor  does  not  have  a  proper  radiometric 
or  wavelength  calibration,  or  when  the  solar  illumination 
intensity  is  unknown.  The  computational  speed  of  the 
atmospheric  correction  method  is  significantly  faster  than  for  the 
first-principles  methods,  making  it  potentially  suitable  for  real¬ 
time  applications.  The  aerosol  optical  depth  retrieval  method, 
unlike  most  prior  methods,  does  not  require  the  presence  of  dark 
pixels.  In  this  paper,  QUAC  is  applied  to  atmospherically 
correction  several  AVIRIS  data  sets.  Comparisons  to  the 
physics-based  FLAASH  code  are  also  presented. 

L  INTRODUCTION 

A  critical  first  step  in  the  analysis  of  visible  through  short- 
wavelength  infrared  (SWIR)  hyperspectral  or  multispectral 
imagery  (HSI  or  MSI)  is  atmospheric  correction  or 
compensation,  in  which  atmospheric  absorption  and  scattering 
effects  are  removed  and  the  data  are  reduced  to  surface 
spectral  reflectance  (see  Fig.  1).  A  number  of  atmospheric- 
correction  methods  and  algorithms  exist,  including  algorithms 
based  on  first-principles  radiation  transport  calculations  [1- 
10],  and  empirical  approaches  such  as  the  Empirical  Line 
Method  (ELM)  [11,12],  which  relies  on  two  or  more  known 
reflectances  in  the  image.  However,  none  of  these  methods 
provides  the  ideal  combination  of  high  accuracy,  high 
computational  speed  and  independence  from  the  need  for  prior 
knowledge. 

This  paper  focuses  on  validation  of  a  new,  semi-empirical 
atmospheric-correction  method,  dubbed  QUAC,  which  also 
enables  retrieval  of  the  wavelength-dependent  optical  depth  of 
aerosol  or  haze  and  molecular  absorbers.  A  more  detailed 


description  of  the  QUAC  approach  was  recently  presented  at 
the  2004  AVIRIS  Workshop  [13]  and  it  is  only  summarized 
here.  It  allows  the  retrieval  of  approximate  reflectance  spectra 
even  when  the  sensor  does  not  have  a  proper  radiometric  or 
wavelength  calibration,  or  when  the  solar  illumination 
intensity  is  unknown,  such  as  when  a  cloud  deck  is  present. 
Computational  speed  is  much  faster  than  for  the  first- 
principles  methods,  making  it  potentially  suitable  for  real-time 
applications.  In  tests  to  date,  QUAC  has  yielded  remarkably 
good  agreement  with  a  state-of-the-art  first-principles 
algorithm.  Like  the  ELM,  QUAC  assumes  a  linear 
relationship  between  spectral  reflectance  and  measured 
radiance,  a  good  approximation  for  most  scenes. 


Figure  1 .  MODTRAN  calculation  of  the  apparent  reflectance  of  a  vegetation 
pixel  as  observed  from  space  with  nadir  viewing,  a  Mid-Latitude  Summer 
model  atmosphere,  and  a  Rural  aerosol  with  visibility=23km. 


The  standard  radiance  equation  may  be  written  as 

W  ,  C{X) 


pj{X)  -  A{X)  + 


\-S{X)<p{X)> 


pj  (A)  + 


l-S(A)<p(A)> 


<p(A)  >, 


(1) 


where  pj  is  the  observed  reflectance  (the  radiance  normalized 
by  the  surface  normal  component  of  the  solar  flux)  for  the  j  ’th 
pixel  at  a  spectral  band  centered  at  wavelength  X  pj°  is  the 
actual  surface  reflectance,  <p>  is  a  spatially  averaged  surface 
reflectance.  A,  B,  C  and  S  are  coefficients  that  describe  the 
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transmission  and  scattering  effects  of  the  atmosphere.  Their 
physical  origin  is  highlighted  in  Figure  2.  The  first 
coefficient,  A,  accounts  for  light  that  never  encounters  the 
surface,  but  is  scattered  and  absorbed  within  the  atmosphere. 
The  second,  B,  accounts  for  the  Sun-surface-sensor  path  direct 
transmittance.  The  third,  C,  accounts  for  diffuse  transmittance 
and  gives  rise  to  the  “adjacency  effect,”  a  spatial  blending 
induced  by  atmospheric  scattering.  The  length  scale  of  the 
adjacency  effect  is  typically  of  order  ~0.5km,  thus  <p>  is 
typically  a  slowly  varying  function  of  position  within  a  large 
image.  S,  the  atmospheric  spherical  albedo,  accounts  for 
enhancement  of  the  ground  illumination  by  atmospheric 
reflection. 


Ptrue  PavB 


Figure  2  .  Radiation-transfer  eontributions  to  the  observed  apparent 
refleetanee,  pobs- 

Eq.  (1)  reduces  to  a  linear  form  under  many  common 
conditions  in  which:  (1)  S<p>  is  small  and  when  either  (2)  the 
diffuse  and  direct  transmittance  terms  can  be  combined  with  a 
single  reflectance  variable,  or  (3)  the  diffuse  term  can  be 
combined  with  the  backscattering  term.  Situation  (1)  occurs 
frequently,  when  the  visibility  is  reasonably  high  or  when  the 
ground  is  dark  in  the  visible  (such  as  with  vegetation,  water  or 
dark  soil).  Situation  (2)  occurs  when  the  pixels  are  very  large, 
several  hundred  meters  in  size.  Situation  (3)  occurs  when  the 
scene  materials  are  fairly  uniformly  interspersed  or  when  the 
image  covers  a  small  geographic  area  (<  ~1  km),  making  <p> 
nearly  constant,  or  when  the  visibility  is  high,  making  the 
diffuse  transmittance  term  small.  When  (1)  and  (3)  apply,  Eq. 
(1)  reduces  to  the  linear  equation 

Pj  (A)  =  A(A)  +  B{X)p°  (A)  +  C(A)  <  p(A)  > .  (2) 

With  the  linear  Eq.  (2),  the  aim  of  atmospheric 
compensation  is  essentially  the  determination  of  an  offset, 
A+C<p>,  and  gain  parameter,  B,  in  order  to  retrieve  the 
surface  reflectance,  Pj°.  Numerous  approaches  to  this  problem 
have  been  developed.  The  ELM  assumes  that  the  radiance 
image  contains  some  pixels  of  known  reflectance.  This 
method  is  not  generally  applicable,  as  in- scene  known 
reflectances  are  often  not  available. 

First-principles  methods  express  the  Eq.  (1)  or  Eq.  (2) 
parameters  in  terms  of  atmospheric  physical  variables,  such  as 
column  water  vapor  and  aerosol  optical  depth  or  visibility. 
For  retrieving  optical  depth,  methods  are  available  that  rely  on 
modeling  the  aerosol  backscatter  over  “dark”  pixels  such  as 
vegetation  and  dark  soil  [14]  or  water  bodies.  However, 
difficulties  in  determining  the  optical  depth  arise  when  there  is 
a  lack  of  suitable  dark  pixels  in  the  scene,  or  when  the  sensor 


is  at  a  low  altitude,  within  the  aerosol  layer,  so  that  the 
backscatter  it  measures  is  a  small  (and  generally  unknown) 
fraction  of  the  total. 

Like  many  first-principles  methods,  QUAC  determines 
the  atmospheric  compensation  parameters  directly  from  the 
information  contained  within  the  scene  (observed  pixel 
spectra),  without  ancillary  information.  However,  unlike  most 
other  methods,  its  aerosol  optical  depth  retrieval  approach 
does  not  require  the  presence  of  dark  pixels.  The  retrieved 
optical  depth  information  can  therefore  be  utilized  to  improve 
the  accuracy  of  methods  that  use  first-principles  modeling.  In 
particular,  it  can  be  used  to  set  the  optical  depth  of  a  model 
aerosol  when  dark  pixels  are  unavailable,  or  to  select  from 
among  alternative  model  aerosols  to  provide  consistency  with 
optical  depths  retrieved  from  a  dark-pixel  method. 

IT  QUAC  ALGORITHM  DESCRIPTION 

The  underlying  assumptions  of  the  approach  are: 

•  There  are  a  number  (~10  or  more)  of  diverse  pixel 
spectra  (diverse  materials)  in  a  scene, 

•  The  spectral  standard  deviation  of  Pj°  for  a  collection 
of  diverse  materials  is  a  nearly  wavelength- 
independent  constant,  and  an  additional,  helpful, 
assumption  is  that, 

•  There  are  sufficiently  dark  pixels  (pj°(?i)~0)  in  a 
scene  to  allow  for  a  good  estimation  of  the  nearly 
spatially  invariant  baseline  contribution, 
Pb=A+C<p>. 

The  first  assumption  is  usually  applicable,  as  it  only 
requires  that  a  handful  of  pixels  out  of  typically  ~10^  to  10^ 
pixels  display  diverse  spectra.  The  most  notable  exception 
would  be  a  scene  over  completely  open  and  deep  water,  in 
which  case  the  material  reflectance  is  well  known  a  priori. 
The  diverse  spectra  can  be  selected  using  any  of  a  number  of 
spectral  diversity  metrics  and  algorithms,  such  as 
endmembers.  The  second  assumption  appears  to  be  generally 
true  based  on  our  empirical  observation,  and  is  likely  related 
to  the  lack  of  spectral  correlation  between  diverse  materials. 
The  third  assumption  is  frequently  applicable,  as  most  scenes 
will  contain  a  number  of  very  dark  pixels  from  such  surfaces 
as  water  bodies,  vegetation  and  cast  shadows.  For  the  atypical 
cases  that  violate  this  assumption,  there  are  alternative 
methods  [13]  for  estimating  a  reasonable  baseline.  The 
atmospheric  correction  step  just  involves  re-arranging  Eq.  (2) 
to  solve  for  Pj°(?L)  given  B  and  the  baseline.  A  key  attribute  of 
QUAC  is  its  applicability  to  any  sensor  viewing  or  solar 
elevation  angle. 

Under  the  above-stated  assumptions,  the  spectral  standard 
deviation  of  Eq.  (2)  for  a  set  of  diverse  pixel  spectral  can  be 
expressed  as 

(jp{X)  =  B{X)ap''  (A) .  (3) 

For  reasons  mentioned  earlier,  A  +  C<p>  can  be  taken  as  a 
constant  in  many,  if  not  most,  cases,  so  it  makes  no 
contribution  to  the  standard  deviation.  In  cases  where  it  varies 


significantly  within  the  scene,  the  image  can  be  divided  into 
smaller  pieces,  as  discussed  below.  Since  op°  is  assumed  to 
be  spectrally  invariant,  then,  to  within  a  normalization  factor 
designated  go,  op  represents  the  correction  factor,  B.  The 
actual  surface  spectral  reflectance  can  be  retrieved  using  the 
extracted  in-scene-determined  compensation  parameters  and 
re-arranging  Eq.  (2)  to  yield 


p0-pb(^) 

go<^Pi^ 


(4) 


where  pb=A+C<p>  is  the  baseline  contribution.  It  is  noted 
that  the  B(X)  is  a  direct  measure  of  the  wavelength-dependent 
aerosol  extinction  and  can  be  used  to  retrieve  the  aerosol 
optical  properties.  The  approach  for  doing  this  was  previously 
described  [13].  Our  focus  here  is  on  the  application  of  QUAC 
to  atmospheric  correction. 


III.  VALIDATION  USING  REAL  AND  SIMULATED 
IMAGERY 


QUAC  was  evaluated  against  a  wide  variety  of  data 
consisting  of  multi-  and  hyperspectral  imagery  for  different 
types  of  scenes,  urban  and  rural,  and  spanning  a  wide  range  of 
atmospheric  conditions.  The  data  are  from  the  airborne 
hyperspectral  AVIRIS  sensor  (224  spectral  channels  from 
-400  to  2500nm,  2-20m  GSD). 

QUAC  was  used  to  perform  atmospheric  correction  and 
aerosol  property  retrieval  on  two  very  different  AVIRIS  data 
collects.  As  depicted  in  Figure  3,  one  corresponds  to  high 
visibility  and  moderate  humidity,  and  the  other  to  low 
visibility  and  high  humidity.  The  NASA  Stennis  data  is 
particularly  useful  because  the  Stennis  site  contains  a  large 
number  of  ground  truth  materials/panels  (these  are  visible  in 
the  lower  left  comer  of  Figure  3). 

The  first  step  in  the  process  is  the  selection  of  diverse 
pixel  spectra.  For  this  analysis,  we  used  the  fast  and 
automated  SMACC  (Sequential  Maximum  Angle  Convex 
Cone)  [16]  endmember  code.  We  used  only  ten  window  region 
wavelengths  in  order  to  further  speed  up  finding  the 
endmembers.  The  results  for  the  NASA  Stennis  scene  are 
displayed  in  Figure  4.  It  is  evident  that  this  set  spans  a  wide 
variety  of  spectral  shapes  and  reflectance  values.  Several 
endmembers  are  quite  dark,  and  the  lowest  reflectance  value 
for  each  channel  defines  the  baseline  spectmm. 

The  next  step  is  to  compute  the  standard  deviation  of  the 
selected  pixels.  Before  this  is  done,  some  refining  of  the 
initial  selection  usually  occurs.  This  involves  weeding  out 
spectra  with  sharp  features,  mainly  vegetation  spectra  that 
display  a  steep  rise  around  700nm  (the  chlorophyll  red  edge). 
Pixels  containing  cirms  clouds,  which  can  be  easily  discerned 
using  established  algorithms,  are  also  rejected.  The  standard 
deviations  for  the  NASA  Stennis  and  North  Carolina  SCAR 
data  are  presented  in  Figure  5.  The  absorption  due  to  the 
940nm  H2O  band  is  clearly  evident,  and  the  much  deeper 
feature  seen  in  the  North  Carolina  data  is  indicative  of  a  much 


higher  humidity  level.  Additional,  weaker  absorption  features, 
such  as  the  840nm  H2O  and  760nm  O2  bands  are  easily 
discernible.  The  general  upper  bounding  envelope  to  these 
curves,  formed  by  spectral  regions  outside  of  the  absorption 
features,  is  a  direct  measure  of  the  aerosol  extinction  for  the  L- 
shaped  path  from  the  Sun  to  the  surface  to  the  sensor  (i.e.,  the 
B  coefficient).  By  inspection,  it  is  quite  obvious  that  the 
Stennis  scene  corresponds  to  a  high  visibility  while  the  North 
Carolina  scene  displays  approximately  an  order  of  magnitude 
more  aerosol  extinction  and  hence  a  much  lower  visibility. 


Figure  3.  RGB  image  of  the  AVIRIS  data  sets  used  for  evaluation  of  QUAC. 
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Figure  4.  The  first  twenty  endmembers  seleeted  by  SMACC  using  the 
apparent  refieetanees  from  the  NASA  Stennis  data. 


Figure  5.  Speetral  standard  deviations  based  on  the  seleeted  endmembers  for 
the  Stennis  and  North  Carolina  AVIRIS  data  sets. 
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From  the  above  information,  baseline  and  standard 
deviation,  the  entire  data  cube  can  be  atmospherically 
corrected.  Sample  results  for  the  Stennis  data  are  shown  in 
Figure  6  and  include  comparisons  to  FLAASH  [7,8]  results 
and  ground  truth  measurements.  In  this  instance,  and  in 
general,  QUAC  compares  well  to  FLAASH.  The 
computational  time  required  for  the  end-to-end  QUAC  process 
for  an  entire  AVIRIS  data  cube  (512x512  pixels  and  224 
spectral  channels  is  -Irnin  on  a  1.6GHz  Pentium  IV  PC.  This 
is  based  on  relatively  slow  IDF  coding  for  the  endmember 
selection  and  atmospheric  correction  steps.  For  comparison, 
FLAASH,  which  run  a  series  of  MODTRAN  calculations, 
requires  ~5-10min  to  perform  the  atmospheric  correction. 
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Figure  6.  Comparison  of  QUAC  atmospherically-corrected  reflectances  to 
those  obtained  from  FLAASH  and  ground  truth  measurements  for  the  Stennis 
data. 

IV.  CONCLUSIONS  AND  FUTURE  DIRECTIONS 

A  new  semi-empirical  algorithm,  QUAC,  for  atmospheric 
correction  and  aerosol  optical  properties  retrieval  for  VNIR- 
SWIR  HSI  and  MSI  sensors  has  been  developed.  Initial 
applications  of  QUAC  to  atmospheric  correction  of  HSI 
AVIRIS  and  MSI  Landsat-7  data  [13]  and  simulated  HSI 
HyMap  data  [13]  show  surprisingly  good  performance,  nearly 
comparable  to  that  of  a  first-principles  physics-based  code. 
Continued  development  and  validation  of  QUAC  is  planned 
using  a  wider  variety  of  HSI  and  MSI  data  sets,  including 
simulated  data,  and  through  field  measurements  involving  full 
characterization  of  the  aerosol  column  concurrent  with 
airborne  and/or  satellite-based  HSI  and  MSI  observations. 
Computational  speed-ups,  automation  and,  eventually,  the 
development  of  an  on-board  data  processing  capability  will 
also  be  explored. 
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